Historical Contingency of Host-Microbiome Response to Perturbation

Infection & Mortality Results

Overview

Setup

(hide)

Click on tabs to display additional information.

Libraries

# Load required libraries for statistical analysis and table creation
library(dplyr)
library(ggplot2)
library(gt)
library(broom)
library(car)
library(emmeans)
library(multcomp)
library(MASS)  # For negative binomial regression
library(rstatix)
library(coin)
library(phyloseq)
library(microViz)
library(tidyverse)

Plotting

# Define treatment order and color palette
treatment_order <- c(
  "A- T- P-",  # Control
  "A- T- P+",  # Parasite
  "A+ T- P-",  # Antibiotics
  "A+ T- P+",  # Antibiotics_Parasite
  "A- T+ P-",  # Temperature
  "A- T+ P+",  # Temperature_Parasite
  "A+ T+ P-",  # Antibiotics_Temperature
  "A+ T+ P+"   # Antibiotics_Temperature_Parasite
)

# Custom color palette matching treatment order
treatment_colors <- c(
  "#1B9E77",  # A- T- P- (Control)
  "#D95F02",  # A- T- P+ (Parasite)
  "#7570B3",  # A+ T- P- (Antibiotics)
  "#E7298A",  # A+ T- P+ (Antibiotics_Parasite)
  "#66A61E",  # A- T+ P- (Temperature)
  "#E6AB02",  # A- T+ P+ (Temperature_Parasite)
  "#A6761D",  # A+ T+ P- (Antibiotics_Temperature)
  "#666666"   # A+ T+ P+ (Antibiotics_Temperature_Parasite)
)

# Create named vector for color scale
treatment_color_scale <- setNames(treatment_colors, treatment_order)

Functions

# Function to extract sample data as dataframe from phyloseq object
samdatAsDataframe <- function(ps) {
  samdat <- phyloseq::sample_data(ps)
  df <- data.frame(samdat, check.names = FALSE, stringsAsFactors = FALSE)
  return(df)
}

# Function to rename variables in phyloseq object
ps_rename <- function(ps, ...) {
  ps <- microViz::ps_get(ps)
  df <- samdatAsDataframe(ps)
  df <- dplyr::rename(.data = df, ...)
  phyloseq::sample_data(ps) <- df
  return(ps)
}

Import Data

ps.unclean <- readRDS('/Users/michaelsieler/Dropbox/Mac (2)/Documents/Sharpton_Lab/Projects_Repository/Rules_of_Life/major-experiment-2023/Data/Robjects/pseq_uncleaned_05052025.rds')


# Create a dataframe for mortality using the "unclean" phyloseq object
# The unclean phyloseq object is prior to removing samples during post-DADA2 processing (e.g., insufficient reads/sample)

data.mortality <-
    ps.unclean %>%

    ## Update Metadata
    ps_rename(Time = Timepoint) %>%
    microViz::ps_mutate(
        Treatment = case_when(
            Antibiotics == 0 & Temperature == 0 & Pathogen == 0 ~ "A- T- P-",
            Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "A- T- P+",
            Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "A+ T- P-",
            Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "A+ T- P+",
            Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "A- T+ P-",
            Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "A- T+ P+",
            Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "A+ T+ P-",
            Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "A+ T+ P+",
            TRUE ~ "Unknown"
        ), .after = "Pathogen"
    ) %>%
    microViz::ps_mutate(Sample = fecal.sample.number, .before = 1) %>%
    microViz::ps_mutate(Sample = gsub("^f", "", Sample)) %>%
    microViz::ps_filter(Treatment != "Unknown") %>%
    microViz::ps_mutate(
        History = case_when(
            Antibiotics + Temperature == 0 ~ 0,
            Antibiotics + Temperature == 1 ~ 1,
            Antibiotics + Temperature == 2 ~ 2,
        ), .after = "Treatment"
    ) %>%
    
    ## Additional metadata updates, factorizing metadata
    microViz::ps_mutate(
    # Create treatment code
        treatment_code = case_when(
          Antibiotics == 0 & Temperature == 0 & Pathogen == 0 ~ "Aneg_Tneg_Pneg",
          Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "Aneg_Tneg_Ppos",
          Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "Apos_Tneg_Pneg",
          Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "Apos_Tneg_Ppos",
          Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "Aneg_Tpos_Pneg",
          Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "Aneg_Tpos_Ppos",
          Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "Apos_Tpos_Pneg",
          Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "Apos_Tpos_Ppos"
        ),
        # Create treatment group factor
        treatment_group = case_when(
          Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "Parasite",
          Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "Antibiotics",
          Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "Antibiotics_Parasite",
          Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "Temperature",
          Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "Temperature_Parasite",
          Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "Antibiotics_Temperature",
          Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "Antibiotics_Temperature_Parasite",
          TRUE ~ "Control"
        ),
        # Convert to factor with appropriate levels
        treatment_group = factor(treatment_group, 
                               levels = c("Control", "Parasite", 
                                          "Antibiotics", "Antibiotics_Parasite",
                                          "Temperature", "Temperature_Parasite",
                                        "Antibiotics_Temperature", "Antibiotics_Temperature_Parasite")
                               ),
        treatment_code = factor(treatment_code, levels = treatment_order),
        # Create time point factor
        time_point = factor(Time, levels = c(0, 14, 18, 25, 29, 60)),
        # Create pathogen status factor
        pathogen_status = factor(ifelse(Pathogen == 1, "Exposed", "Unexposed"),
                               levels = c("Unexposed", "Exposed")),
        # Create sex factor
        sex = factor(Sex, levels = c("M", "F"))
        ) %>%
    microViz::ps_mutate(Treatment = factor(Treatment, levels = treatment_order)) %>%
      microViz::ps_mutate(Exp_Type = case_when(
          Treatment %in% c("A- T- P-", "A- T- P+")  ~ "No prior stressor(s)",
          Treatment %in% c("A+ T- P-", "A+ T- P+")  ~ "Antibiotics",
          Treatment %in% c("A- T+ P-", "A- T+ P+") ~ "Temperature",
          Treatment %in% c("A+ T+ P-", "A+ T+ P+") ~ "Combined",
      )) %>%
      microViz::ps_mutate(Exp_Type = factor(Exp_Type, levels = c("No prior stressor(s)", "Antibiotics", "Temperature", "Combined"))) %>%
    microViz::samdat_tbl() 


data.infection <- data.mortality %>%
    dplyr::filter(Pathogen == 1) %>%
    dplyr::group_by(Treatment) %>%
    dplyr::summarise(
        n_fish = n(),
        n_infected = sum(Total.Worm.Count > 0),
        mean_worm_burden = mean(Total.Worm.Count, na.rm = TRUE),
        percent_infected = round((n_infected / n_fish) * 100, 1),
        .groups = "drop"
    ) %>%
    dplyr::mutate(
        # Convert Treatment to factor with specific order
        Treatment = factor(Treatment, levels = treatment_order)
    )

data.infection_tank <- data.mortality %>%
    dplyr::filter(Pathogen == 1) %>%
    dplyr::group_by(Treatment, Tank.ID) %>%
    dplyr::summarise(
        n_fish = n(),
        n_infected = sum(Total.Worm.Count > 0),
        mean_worm_burden = mean(Total.Worm.Count, na.rm = TRUE),
        percent_infected = round((n_infected / n_fish) * 100, 1),
        .groups = "drop"
    ) %>%
    dplyr::mutate(
        # Convert Treatment to factor with specific order
        Treatment = factor(Treatment, levels = treatment_order)
    )

# data.mortality
# data.infection
# data.infection_tank

00) All

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality

All

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

By Tank

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

Statistical Analysis - All Treatments Mortality

Code

# Prepare data for statistical analysis - percent mortality rates
mortality_stats_data <- data.mortality %>%
  dplyr::group_by(Treatment) %>%
  dplyr::count() %>%
  dplyr::mutate(
    total_per_group = 90,
    dead = total_per_group - n,
    percent_mortality = round((dead / total_per_group) * 100, 1)
  )

# Chi-square test for independence between treatment and mortality
set.seed(42)
chi_square_test <- chisq.test(
  matrix(c(mortality_stats_data$dead, mortality_stats_data$n), 
         ncol = 2, byrow = FALSE)
)

# Create mortality rate summary table for display
mortality_rate_table <- mortality_stats_data %>%
  dplyr::select(Treatment, dead, n, percent_mortality) %>%
  dplyr::rename(Survived = n) %>%
  dplyr::mutate(
    Total = dead + Survived,
    percent_survived = round((Survived / Total) * 100, 1)
  ) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Mortality Rate Summary - All Treatments",
    subtitle = "Percent mortality and survival rates by treatment"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(dead, Survived, Total),
    decimals = 0
  ) %>%
  gt::fmt_number(
    columns = c(percent_mortality, percent_survived),
    decimals = 1
  )

# Chi-square test results
test_stat_table <- data.frame(
  Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
  Value = c(
    round(chi_square_test$statistic, 4),
    chi_square_test$parameter,
    round(chi_square_test$p.value, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Chi-Square Test Results - All Treatments Mortality",
    subtitle = "Test for independence between treatment and mortality"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Post-hoc pairwise comparisons using chi-square tests
set.seed(42)
pairwise_results <- list()
treatments <- mortality_stats_data$Treatment

for(i in 1:(length(treatments)-1)) {
  for(j in (i+1):length(treatments)) {
    trt1 <- treatments[i]
    trt2 <- treatments[j]
    
    # Get data for these two treatments
    data1 <- mortality_stats_data %>% dplyr::filter(Treatment == trt1)
    data2 <- mortality_stats_data %>% dplyr::filter(Treatment == trt2)
    
    # Create 2x2 contingency table
    cont_table <- matrix(
      c(data1$dead, data1$n, data2$dead, data2$n),
      nrow = 2, ncol = 2,
      dimnames = list(
        c("Dead", "Survived"),
        c(trt1, trt2)
      )
    )
    
    # Perform chi-square test
    test_result <- chisq.test(cont_table)
    
    pairwise_results[[paste(trt1, "vs", trt2)]] <- data.frame(
      Comparison = paste(trt1, "vs", trt2),
      Chi_square = round(test_result$statistic, 4),
      P_value = round(test_result$p.value, 6),
      Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
    )
  }
}

# Combine pairwise results
pairwise_summary <- do.call(rbind, pairwise_results) %>%
    dplyr::arrange(desc(Significant), P_value) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pairwise Chi-Square Tests - All Treatments Mortality",
    subtitle = "Post-hoc comparisons between treatment pairs"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )

# Show only significant pairwise comparisons
pairwise_summary__SigOnly <- do.call(rbind, pairwise_results) %>%
    dplyr::arrange(desc(Significant), P_value) %>%
    dplyr::filter(Significant == "Yes") %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Significant Pairwise Chi-Square Tests - All Treatments Mortality",
    subtitle = "Post-hoc comparisons between treatment pairs (significant only)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )

Table

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality Rates

Mortality Rate Summary - All Treatments
Percent mortality and survival rates by treatment
dead	Survived	percent_mortality	Total	percent_survived
A- T- P-
9	81	10.0	90	90.0
A- T- P+
8	82	8.9	90	91.1
A+ T- P-
14	76	15.6	90	84.4
A+ T- P+
10	80	11.1	90	88.9
A- T+ P-
19	71	21.1	90	78.9
A- T+ P+
13	77	14.4	90	85.6
A+ T+ P-
30	60	33.3	90	66.7
A+ T+ P+
20	70	22.2	90	77.8

Chi-Square

Chi-Square Test Results - All Treatments Mortality
Test for independence between treatment and mortality
Statistic	Value
Chi-square statistic	29.797800
Degrees of freedom	7.000000
P-value	0.000103

Chi-Square Pairwise

Pairwise Chi-Square Tests - All Treatments Mortality
Post-hoc comparisons between treatment pairs
Comparison	Chi_square	P_value	Significant
A- T- P+ vs A+ T+ P-	14.7109	0.000125	Yes
A- T- P- vs A+ T+ P-	13.0933	0.000296	Yes
A+ T- P+ vs A+ T+ P-	11.6036	0.000658	Yes
A- T+ P+ vs A+ T+ P-	7.8221	0.005161	Yes
A+ T- P- vs A+ T+ P-	6.7680	0.009280	Yes
A- T- P+ vs A+ T+ P+	5.1175	0.023686	Yes
A- T- P+ vs A- T+ P-	4.3573	0.036851	Yes
A- T- P- vs A+ T+ P+	4.1105	0.042617	Yes
A- T- P- vs A- T+ P-	3.4258	0.064187	No
A+ T- P+ vs A+ T+ P+	3.2400	0.071861	No
A- T+ P- vs A+ T+ P-	2.8042	0.094019	No
A+ T- P+ vs A- T+ P-	2.6307	0.104813	No
A+ T+ P- vs A+ T+ P+	2.2431	0.134214	No
A- T+ P+ vs A+ T+ P+	1.3358	0.247775	No
A- T- P+ vs A+ T- P-	1.2946	0.255204	No
A- T+ P- vs A- T+ P+	0.9502	0.329676	No
A+ T- P- vs A+ T+ P+	0.9065	0.341038	No
A- T- P+ vs A- T+ P+	0.8625	0.353031	No
A- T- P- vs A+ T- P-	0.7976	0.371823	No
A+ T- P- vs A- T+ P-	0.5937	0.440995	No
A- T- P- vs A- T+ P+	0.4661	0.494809	No
A+ T- P- vs A+ T- P+	0.4327	0.510671	No
A+ T- P+ vs A- T+ P+	0.1994	0.655213	No
A- T- P+ vs A+ T- P+	0.0617	0.803785	No
A- T- P- vs A- T- P+	0.0000	1.000000	No
A- T- P- vs A+ T- P+	0.0000	1.000000	No
A+ T- P- vs A- T+ P+	0.0000	1.000000	No
A- T+ P- vs A+ T+ P+	0.0000	1.000000	No

Significant Pairwise Chi-Square Tests - All Treatments Mortality
Post-hoc comparisons between treatment pairs (significant only)
Comparison	Chi_square	P_value	Significant
A- T- P+ vs A+ T+ P-	14.7109	0.000125	Yes
A- T- P- vs A+ T+ P-	13.0933	0.000296	Yes
A+ T- P+ vs A+ T+ P-	11.6036	0.000658	Yes
A- T+ P+ vs A+ T+ P-	7.8221	0.005161	Yes
A+ T- P- vs A+ T+ P-	6.7680	0.009280	Yes
A- T- P+ vs A+ T+ P+	5.1175	0.023686	Yes
A- T- P+ vs A- T+ P-	4.3573	0.036851	Yes
A- T- P- vs A+ T+ P+	4.1105	0.042617	Yes

01) What is the host-microbiome response to parasite exposure?

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality

All

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

By Tank

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

Statistical Analysis - Parasite Exposure Mortality

Code

# Prepare data for statistical analysis - parasite exposure comparison (percent mortality rates)
parasite_mortality_data <- data.mortality %>%
  dplyr::filter(Treatment %in% c("A- T- P-", "A- T- P+")) %>%
  dplyr::group_by(Treatment) %>%
  dplyr::count() %>%
  dplyr::mutate(
    total_per_group = 90,
    dead = total_per_group - n,
    percent_mortality = round((dead / total_per_group) * 100, 1)
  )

# Fisher's exact test for 2x2 contingency table (more appropriate for small samples)
set.seed(42)
fisher_test <- fisher.test(
  matrix(c(parasite_mortality_data$dead, parasite_mortality_data$n), 
         ncol = 2, byrow = FALSE)
)

# Chi-square test as alternative
chi_square_parasite <- chisq.test(
  matrix(c(parasite_mortality_data$dead, parasite_mortality_data$n), 
         ncol = 2, byrow = FALSE)
)

# Create mortality rate summary table for display
parasite_mortality_rate_table <- parasite_mortality_data %>%
  dplyr::select(Treatment, dead, n, percent_mortality) %>%
  dplyr::rename(Survived = n) %>%
  dplyr::mutate(
    Total = dead + Survived,
    percent_survived = round((Survived / Total) * 100, 1)
  ) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Mortality Rate Summary - Parasite Exposure",
    subtitle = "Percent mortality and survival rates: A-T-P- vs A-T-P+"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(dead, Survived, Total),
    decimals = 0
  ) %>%
  gt::fmt_number(
    columns = c(percent_mortality, percent_survived),
    decimals = 1
  )

# Fisher's exact test results
fisher_results <- data.frame(
  Statistic = c("Fisher's exact test p-value", "Odds ratio", "95% CI lower", "95% CI upper"),
  Value = c(
    round(fisher_test$p.value, 6),
    round(fisher_test$estimate, 4),
    round(fisher_test$conf.int[1], 4),
    round(fisher_test$conf.int[2], 4)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Fisher's Exact Test Results - Parasite Exposure Mortality",
    subtitle = "Comparison of A-T-P- vs A-T-P+ mortality rates"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Chi-square test results for comparison
chi_square_parasite_results <- data.frame(
  Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
  Value = c(
    round(chi_square_parasite$statistic, 4),
    chi_square_parasite$parameter,
    round(chi_square_parasite$p.value, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Chi-Square Test Results - Parasite Exposure Mortality",
    subtitle = "Alternative test for A-T-P- vs A-T-P+ mortality rates"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Additional analysis: Compare mortality rates between parasite-exposed and unexposed
# Calculate difference in mortality rates
parasite_mortality_diff <- parasite_mortality_data %>%
  dplyr::arrange(Treatment) %>%
  dplyr::mutate(
    Pathogen_Status = dplyr::case_when(
      stringr::str_detect(Treatment, "P\\+") ~ "Exposed",
      stringr::str_detect(Treatment, "P\\-") ~ "Unexposed",
      TRUE ~ "Unknown"
    )
  ) %>%
  dplyr::select(Pathogen_Status, percent_mortality)

## Adding missing grouping variables: `Treatment`

# Calculate the difference in mortality rates
mortality_rate_diff <- parasite_mortality_diff$percent_mortality[2] - parasite_mortality_diff$percent_mortality[1]

# Create summary of mortality rate difference
parasite_mortality_diff_summary <- data.frame(
  Comparison = "A-T-P+ vs A-T-P-",
  Mortality_Rate_Difference = round(mortality_rate_diff, 2),
  Unexposed_Rate = round(parasite_mortality_diff$percent_mortality[1], 2),
  Exposed_Rate = round(parasite_mortality_diff$percent_mortality[2], 2),
  Percent_Increase = round((mortality_rate_diff / parasite_mortality_diff$percent_mortality[1]) * 100, 1)
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Mortality Rate Difference - Parasite Exposure",
    subtitle = "Comparison of mortality rates between parasite-exposed and unexposed groups"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(Mortality_Rate_Difference, Unexposed_Rate, Exposed_Rate, Percent_Increase),
    decimals = 2
  )

Table

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality Rates

Mortality Rate Summary - Parasite Exposure
Percent mortality and survival rates: A-T-P- vs A-T-P+
dead	Survived	percent_mortality	Total	percent_survived
A- T- P-
9	81	10.0	90	90.0
A- T- P+
8	82	8.9	90	91.1

Mortality Rate Difference

Mortality Rate Difference - Parasite Exposure
Comparison of mortality rates between parasite-exposed and unexposed groups
Comparison	Mortality_Rate_Difference	Unexposed_Rate	Exposed_Rate	Percent_Increase
A-T-P+ vs A-T-P-	−1.10	10.00	8.90	−11.00

Fisher’s Exact

Fisher's Exact Test Results - Parasite Exposure Mortality
Comparison of A-T-P- vs A-T-P+ mortality rates
Statistic	Value
Fisher's exact test p-value	1.0000
Odds ratio	1.1381
95% CI lower	0.3692
95% CI upper	3.5739

Chi-Square

Chi-Square Test Results - Parasite Exposure Mortality
Alternative test for A-T-P- vs A-T-P+ mortality rates
Statistic	Value
Chi-square statistic	0
Degrees of freedom	1
P-value	1

02) Is the host-microbiome response to parasite exposure historically contingent?

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality

All

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

By Tank

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

Statistical Analysis - Historical Contingency Mortality

Code

# Prepare data for statistical analysis - historical contingency mortality
historical_mortality_data <- data.mortality %>%
  dplyr::filter(Treatment %in% c("A- T- P+", "A+ T- P+", "A- T+ P+", "A+ T+ P+")) %>%
  dplyr::group_by(Treatment) %>%
  dplyr::count() %>%
  dplyr::mutate(
    total_per_group = 90,
    dead = total_per_group - n,
    percent_mortality = round((dead / total_per_group) * 100, 1)
  )

# Chi-square test for independence between treatment and mortality
set.seed(42)
historical_chi_square_test <- chisq.test(
  matrix(c(historical_mortality_data$dead, historical_mortality_data$n), 
         ncol = 2, byrow = FALSE)
)

# Create mortality rate summary table for display
historical_mortality_rate_table <- historical_mortality_data %>%
  dplyr::select(Treatment, dead, n, percent_mortality) %>%
  dplyr::rename(Survived = n) %>%
  dplyr::mutate(
    Total = dead + Survived,
    percent_survived = round((Survived / Total) * 100, 1)
  ) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Mortality Rate Summary - Historical Contingency",
    subtitle = "Percent mortality and survival rates by treatment (parasite-exposed only)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(dead, Survived, Total),
    decimals = 0
  ) %>%
  gt::fmt_number(
    columns = c(percent_mortality, percent_survived),
    decimals = 1
  )

# Chi-square test results
historical_chi_square_results <- data.frame(
  Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
  Value = c(
    round(historical_chi_square_test$statistic, 4),
    historical_chi_square_test$parameter,
    round(historical_chi_square_test$p.value, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Chi-Square Test Results - Historical Contingency Mortality",
    subtitle = "Test for independence between treatment and mortality (parasite-exposed only)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Post-hoc pairwise comparisons using chi-square tests
set.seed(42)
historical_pairwise_results <- list()
historical_treatments <- historical_mortality_data$Treatment

for(i in 1:(length(historical_treatments)-1)) {
  for(j in (i+1):length(historical_treatments)) {
    trt1 <- historical_treatments[i]
    trt2 <- historical_treatments[j]
    
    # Get data for these two treatments
    data1 <- historical_mortality_data %>% dplyr::filter(Treatment == trt1)
    data2 <- historical_mortality_data %>% dplyr::filter(Treatment == trt2)
    
    # Create 2x2 contingency table
    cont_table <- matrix(
      c(data1$dead, data1$n, data2$dead, data2$n),
      nrow = 2, ncol = 2,
      dimnames = list(
        c("Dead", "Survived"),
        c(trt1, trt2)
      )
    )
    
    # Perform chi-square test
    test_result <- chisq.test(cont_table)
    
    historical_pairwise_results[[paste(trt1, "vs", trt2)]] <- data.frame(
      Comparison = paste(trt1, "vs", trt2),
      Chi_square = round(test_result$statistic, 4),
      P_value = round(test_result$p.value, 6),
      Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
    )
  }
}

# Combine pairwise results
historical_pairwise_summary <- do.call(rbind, historical_pairwise_results) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pairwise Chi-Square Tests - Historical Contingency Mortality",
    subtitle = "Post-hoc comparisons between treatment pairs (parasite-exposed only)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )


data.mortality %>%
    dplyr::filter(Treatment %in% c("A- T- P+", "A+ T- P+", "A- T+ P+", "A+ T+ P+")) %>%
  dplyr::group_by(Treatment) %>%
  dplyr::count() %>%
  dplyr::mutate(
    total_per_group = 90,
    dead = total_per_group - n,
    percent_mortality = round((dead / total_per_group) * 100, 1)
  )

## # A tibble: 4 × 5
## # Groups:   Treatment [4]
##   Treatment     n total_per_group  dead percent_mortality
##   <fct>     <int>           <dbl> <dbl>             <dbl>
## 1 A- T- P+     82              90     8               8.9
## 2 A+ T- P+     80              90    10              11.1
## 3 A- T+ P+     77              90    13              14.4
## 4 A+ T+ P+     70              90    20              22.2

Table

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality Rates

Mortality Rate Summary - Historical Contingency
Percent mortality and survival rates by treatment (parasite-exposed only)
dead	Survived	percent_mortality	Total	percent_survived
A- T- P+
8	82	8.9	90	91.1
A+ T- P+
10	80	11.1	90	88.9
A- T+ P+
13	77	14.4	90	85.6
A+ T+ P+
20	70	22.2	90	77.8

Chi-Square

Chi-Square Test Results - Historical Contingency Mortality
Test for independence between treatment and mortality (parasite-exposed only)
Statistic	Value
Chi-square statistic	7.561400
Degrees of freedom	3.000000
P-value	0.056002

Chi-Square Pairwise

Pairwise Chi-Square Tests - Historical Contingency Mortality
Post-hoc comparisons between treatment pairs (parasite-exposed only)
Comparison	Chi_square	P_value	Significant
A- T- P+ vs A+ T- P+	0.0617	0.803785	No
A- T- P+ vs A- T+ P+	0.8625	0.353031	No
A- T- P+ vs A+ T+ P+	5.1175	0.023686	Yes
A+ T- P+ vs A- T+ P+	0.1994	0.655213	No
A+ T- P+ vs A+ T+ P+	3.2400	0.071861	No
A- T+ P+ vs A+ T+ P+	1.3358	0.247775	No

Infection

All

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

By Tank

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

Statistical Analysis - Infection Prevalence and Worm Burden

Code

# Statistical analysis for infection prevalence and worm burden
set.seed(42)

# 1. INFECTION PREVALENCE ANALYSIS
# Chi-square test for infection prevalence across treatments
infection_chi_square <- chisq.test(
  data.infection %>%
    dplyr::select(n_infected, n_fish) %>%
    dplyr::mutate(n_uninfected = n_fish - n_infected) %>%
    dplyr::select(n_infected, n_uninfected) %>%
    as.matrix()
)

# Create infection prevalence contingency table
infection_contingency_table <- data.infection %>%
  dplyr::mutate(
    n_uninfected = n_fish - n_infected,
    percent_uninfected = round((n_uninfected / n_fish) * 100, 1)
  ) %>%
  dplyr::select(Treatment, n_infected, n_uninfected, n_fish, percent_infected, percent_uninfected) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Infection Prevalence Contingency Table",
    subtitle = "Number of infected vs uninfected fish by treatment"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(n_infected, n_uninfected, n_fish),
    decimals = 0
  )

# Chi-square test results for infection prevalence
infection_chi_results <- data.frame(
  Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
  Value = c(
    round(infection_chi_square$statistic, 4),
    infection_chi_square$parameter,
    round(infection_chi_square$p.value, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Chi-Square Test Results - Infection Prevalence",
    subtitle = "Test for independence between treatment and infection status"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Post-hoc pairwise chi-square tests for infection prevalence
set.seed(42)
infection_pairwise_results <- list()
infection_treatments <- data.infection$Treatment

for(i in 1:(length(infection_treatments)-1)) {
  for(j in (i+1):length(infection_treatments)) {
    trt1 <- infection_treatments[i]
    trt2 <- infection_treatments[j]
    
    # Get data for these two treatments
    data1 <- data.infection %>% dplyr::filter(Treatment == trt1)
    data2 <- data.infection %>% dplyr::filter(Treatment == trt2)
    
    # Create 2x2 contingency table
    cont_table <- matrix(
      c(data1$n_infected, data1$n_fish - data1$n_infected, 
        data2$n_infected, data2$n_fish - data2$n_infected),
      nrow = 2, ncol = 2,
      dimnames = list(
        c("Infected", "Uninfected"),
        c(trt1, trt2)
      )
    )
    
    # Perform chi-square test
    test_result <- chisq.test(cont_table)
    
    infection_pairwise_results[[paste(trt1, "vs", trt2)]] <- data.frame(
      Comparison = paste(trt1, "vs", trt2),
      Chi_square = round(test_result$statistic, 4),
      P_value = round(test_result$p.value, 6),
      Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
    )
  }
}

# Combine pairwise results for infection prevalence
infection_pairwise_summary <- do.call(rbind, infection_pairwise_results) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pairwise Chi-Square Tests - Infection Prevalence",
    subtitle = "Post-hoc comparisons between treatment pairs"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )

# 2. WORM BURDEN ANALYSIS (TOTAL WORMS)
# Prepare individual-level data for negative binomial regression
set.seed(42)
individual_worm_data <- data.mortality %>%
  dplyr::filter(Pathogen == 1) %>%
  dplyr::select(Treatment, Total.Worm.Count) %>%
  dplyr::filter(!is.na(Total.Worm.Count)) %>%
  dplyr::mutate(
    Treatment = factor(Treatment, levels = c("A- T- P+", "A+ T- P+", "A- T+ P+", "A+ T+ P+"))
  )

# Fit negative binomial model for total worm counts
# Using MASS::glm.nb for negative binomial regression
library(MASS)
nb_model <- MASS::glm.nb(Total.Worm.Count ~ Treatment, data = individual_worm_data)

# Get model summary
nb_summary <- summary(nb_model)

nb_model.anova <- Anova(nb_model, type = 2)

# Extract coefficients and significance
nb_coefficients <- data.frame(
  Treatment = c("Intercept", levels(individual_worm_data$Treatment)[-1]),
  Estimate = round(nb_summary$coefficients[, 1], 4),
  Std_Error = round(nb_summary$coefficients[, 2], 4),
  z_value = round(nb_summary$coefficients[, 3], 4),
  p_value = round(nb_summary$coefficients[, 4], 6)
) %>%
  dplyr::mutate(
    p_signif = ifelse(p_value < 0.001, "***", 
                     ifelse(p_value < 0.01, "**",
                            ifelse(p_value < 0.05, "*", "ns")))
  )

# Create worm burden summary table
worm_burden_data <- data.mortality %>%
  dplyr::filter(Pathogen == 1) %>%
  dplyr::group_by(Treatment) %>%
  dplyr::summarise(
    total_worms = sum(Total.Worm.Count, na.rm = TRUE),
    mean_worms = mean(Total.Worm.Count, na.rm = TRUE),
    median_worms = median(Total.Worm.Count, na.rm = TRUE),
    sd_worms = sd(Total.Worm.Count, na.rm = TRUE),
    n_fish = n(),
    .groups = "drop"
  ) %>%
  dplyr::mutate(
    se_worms = sd_worms / sqrt(n_fish)
  )

worm_burden_summary <- worm_burden_data %>%
  dplyr::select(Treatment, total_worms, mean_worms, median_worms, sd_worms, se_worms, n_fish) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Worm Burden Summary by Treatment",
    subtitle = "Total, mean, median, standard deviation, and standard error of worm counts per treatment"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(total_worms, n_fish),
    decimals = 0
  ) %>%
  gt::fmt_number(
    columns = c(mean_worms, median_worms, sd_worms, se_worms),
    decimals = 2
  )

# Negative binomial model results
nb_model_results <- nb_coefficients %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Negative Binomial Regression Results - Worm Burden",
    subtitle = "Model: glm.nb(worm counts ~ Treatment)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(Estimate, Std_Error, z_value, p_value),
    decimals = 3
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = p_signif,
      rows = p_signif != "ns"
    )
  )


nb_model.anova %>%
    broom::tidy() %>%
  dplyr::mutate(
    p_signif = ifelse(p.value < 0.001, "***", 
                     ifelse(p.value < 0.01, "**",
                            ifelse(p.value < 0.05, "*", "ns")))
  ) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Negative Binomial Regression ANOVA Results - Worm Burden",
    subtitle = "Model: Anova(glm.nb(worm counts ~ Treatment))"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(p.value),
    decimals = 3
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = p_signif,
      rows = p_signif != "ns"
    )
  )

Negative Binomial Regression ANOVA Results - Worm Burden
Model: Anova(glm.nb(worm counts ~ Treatment))
term	statistic	df	p.value	p_signif
Treatment	9.044622	3	0.029	*

# Model fit statistics
nb_fit_stats <- data.frame(
  Statistic = c("Log-likelihood", "AIC", "Theta", "SE Theta", "Deviance", "Residual Deviance"),
  Value = c(
    round(logLik(nb_model), 2),
    round(AIC(nb_model), 2),
    round(nb_model$theta, 4),
    round(nb_model$SE.theta, 4),
    round(nb_summary$deviance, 2),
    round(sum(residuals(nb_model, type = "deviance")^2), 2)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Negative Binomial Model Fit Statistics",
    subtitle = "Model diagnostics and fit measures"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Post-hoc pairwise comparisons for worm burden using emmeans
set.seed(42)
library(emmeans)

# Get estimated marginal means and pairwise comparisons
nb_emmeans <- emmeans::emmeans(nb_model, ~ Treatment)
nb_pairwise <- pairs(nb_emmeans, adjust = "bonferroni")

# Convert to data frame for display
nb_pairwise_df <- as.data.frame(nb_pairwise) %>%
  dplyr::mutate(
    contrast = as.character(contrast),
    estimate = round(estimate, 4),
    SE = round(SE, 4),
    z.ratio = round(z.ratio, 4),
    p.value = round(p.value, 6),
    p_signif = ifelse(p.value < 0.001, "***", 
                     ifelse(p.value < 0.01, "**",
                            ifelse(p.value < 0.05, "*", "ns")))
  ) %>%
  dplyr::select(contrast, estimate, SE, z.ratio, p.value, p_signif)

# Display pairwise comparisons for worm burden
nb_pairwise_summary <- nb_pairwise_df %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pairwise Comparisons - Worm Burden (Negative Binomial)",
    subtitle = "Bonferroni-adjusted pairwise comparisons between treatments"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(estimate, SE, z.ratio, p.value),
    decimals = 3
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = p_signif,
      rows = p_signif != "ns"
    )
  )

# Model diagnostics for negative binomial
# Check dispersion parameter and model adequacy
nb_diagnostics <- data.frame(
  Statistic = c("Theta (dispersion parameter)", "SE Theta", "Deviance", "Degrees of Freedom", "Deviance/DF", "Model adequacy"),
  Value = c(
    round(nb_model$theta, 4),
    round(nb_model$SE.theta, 4),
    round(nb_summary$deviance, 2),
    nb_summary$df.residual,
    round(nb_summary$deviance / nb_summary$df.residual, 2),
    ifelse(nb_summary$deviance / nb_summary$df.residual < 2, "Good fit", "Consider model diagnostics")
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Negative Binomial Model Diagnostics",
    subtitle = "Checking model adequacy and dispersion"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Compare with Poisson model for reference
set.seed(42)
poisson_model <- glm(Total.Worm.Count ~ Treatment, data = individual_worm_data, family = "poisson")
poisson_summary <- summary(poisson_model)

# Model comparison
model_comparison <- data.frame(
  Model = c("Poisson", "Negative Binomial"),
  Log_Likelihood = c(round(logLik(poisson_model), 2), round(logLik(nb_model), 2)),
  AIC = c(round(AIC(poisson_model), 2), round(AIC(nb_model), 2)),
  Deviance = c(round(poisson_summary$deviance, 2), round(nb_summary$deviance, 2)),
  DF = c(poisson_summary$df.residual, nb_summary$df.residual),
  Deviance_DF_Ratio = c(
    round(poisson_summary$deviance / poisson_summary$df.residual, 2),
    round(nb_summary$deviance / nb_summary$df.residual, 2)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Model Comparison: Poisson vs Negative Binomial",
    subtitle = "Comparing model fit and adequacy"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(Log_Likelihood, AIC, Deviance, Deviance_DF_Ratio),
    decimals = 2
  ) %>%
  gt::fmt_number(
    columns = DF,
    decimals = 0
  )

# 3. NORMALIZED WORM BURDEN ANALYSIS (WORMS PER INFECTED FISH)
# Calculate worms per infected fish to normalize for infection prevalence
# NOTE: This is ratio data, not count data, so we need different statistical approaches
set.seed(42)

# Create normalized worm burden summary table
normalized_worm_summary <- normalized_worm_data %>%
  dplyr::select(Treatment, total_worms, n_infected, worms_per_infected, n_fish) %>%
  dplyr::mutate(
    infection_rate = round((n_infected / n_fish) * 100, 1)
  ) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Normalized Worm Burden Summary by Treatment",
    subtitle = "Worms per infected fish (normalized for infection prevalence)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(total_worms, n_infected, n_fish),
    decimals = 0
  ) %>%
  gt::fmt_number(
    columns = c(worms_per_infected, infection_rate),
    decimals = 2
  )

# Create tank_normalized_for_test dataset for statistical analysis
# Using ANOVA on worms per infected fish (excluding tanks with 0 infected fish)
set.seed(42)
tank_normalized_for_test <- tank_normalized_data %>%
  dplyr::filter(n_infected > 0) %>%
  dplyr::mutate(
    Treatment = factor(Treatment, levels = treatment_order)
  )

# APPROPRIATE STATISTICAL ANALYSIS FOR RATIO DATA
# Method 1: One-way ANOVA on tank-level ratios (if assumptions are met)
set.seed(42)
# Use tank-level data for statistical testing
normalized_anova <- aov(worms_per_infected ~ Treatment, data = tank_normalized_for_test)
normalized_anova_summary <- summary(normalized_anova)

# Extract ANOVA results using broom package for cleaner extraction
library(broom)
normalized_anova_tidy <- tidy(normalized_anova)

# Calculate F-statistic and p-value manually from ANOVA summary
f_stat <- normalized_anova_summary[[1]]$"Mean Sq"[1] / normalized_anova_summary[[1]]$"Mean Sq"[2]
p_val <- pf(f_stat, normalized_anova_summary[[1]]$"Df"[1], normalized_anova_summary[[1]]$"Df"[2], lower.tail = FALSE)

normalized_anova_results <- data.frame(
  Statistic = c("F-value", "Degrees of freedom (Treatment)", "Degrees of freedom (Residuals)", "P-value"),
  Value = c(
    round(f_stat, 4),
    normalized_anova_summary[[1]]$"Df"[1],
    normalized_anova_summary[[1]]$"Df"[2],
    round(p_val, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "ANOVA Results - Normalized Worm Burden",
    subtitle = "One-way ANOVA: worms per infected fish ~ Treatment (tank-level data)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Post-hoc Tukey test for normalized worm burden
set.seed(42)
normalized_tukey <- TukeyHSD(normalized_anova)

# Convert Tukey results to data frame
normalized_tukey_df <- as.data.frame(normalized_tukey$Treatment) %>%
  dplyr::mutate(
    comparison = rownames(.),
    diff = round(diff, 4),
    lwr = round(lwr, 4),
    upr = round(upr, 4),
    p_adj = round(`p adj`, 6),
    p_signif = ifelse(p_adj < 0.001, "***",
                     ifelse(p_adj < 0.01, "**",
                            ifelse(p_adj < 0.05, "*", "ns")))
  ) %>%
  dplyr::select(comparison, diff, lwr, upr, p_adj, p_signif)

# Display Tukey test results
normalized_tukey_summary <- normalized_tukey_df %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Tukey HSD Test Results - Normalized Worm Burden",
    subtitle = "Post-hoc pairwise comparisons for worms per infected fish (tank-level data)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(diff, lwr, upr, p_adj),
    decimals = 3
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = p_signif,
      rows = p_signif != "ns"
    )
  )

# Method 2: Kruskal-Wallis test (non-parametric alternative)
set.seed(42)
normalized_kruskal <- kruskal.test(worms_per_infected ~ Treatment, data = tank_normalized_for_test)

# Kruskal-Wallis results
normalized_kruskal_results <- data.frame(
  Statistic = c("Kruskal-Wallis chi-squared", "Degrees of freedom", "P-value"),
  Value = c(
    round(normalized_kruskal$statistic, 4),
    normalized_kruskal$parameter,
    round(normalized_kruskal$p.value, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Kruskal-Wallis Test Results - Normalized Worm Burden",
    subtitle = "Non-parametric test for differences in worms per infected fish across treatments (tank-level data)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Method 3: Permutation test (most appropriate for ratio data with small sample sizes)
set.seed(42)
library(coin)

# Create permutation test for normalized worm burden using tank-level data
normalized_permutation <- coin::oneway_test(
  worms_per_infected ~ Treatment,
  data = tank_normalized_for_test,
  distribution = approximate(nresample = 10000)
)

# Permutation test results
normalized_permutation_results <- data.frame(
  Statistic = c("Z statistic", "P-value (approximate)", "Number of permutations"),
  Value = c(
    round(statistic(normalized_permutation), 4),
    round(pvalue(normalized_permutation), 6),
    attr(pvalue(normalized_permutation), "nresample")
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Permutation Test Results - Normalized Worm Burden",
    subtitle = "Non-parametric permutation test for ratio data (tank-level data)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Method 4: Bootstrap confidence intervals for ratios
set.seed(42)
library(boot)

## 
## Attaching package: 'boot'

## The following object is masked from 'package:survival':
## 
##     aml

## The following object is masked from 'package:car':
## 
##     logit

# Function to calculate mean ratio for bootstrap
ratio_mean <- function(data, indices) {
  d <- data[indices, ]
  return(mean(d$worms_per_infected, na.rm = TRUE))
}

# Bootstrap analysis for each treatment using tank-level data
# We need multiple observations per treatment for bootstrap to work
bootstrap_results <- list()
for(trt in unique(tank_normalized_data$Treatment)) {
  trt_data <- tank_normalized_data %>% dplyr::filter(Treatment == trt, n_infected > 0)

  # Only proceed if we have multiple observations
  if(nrow(trt_data) > 1) {
    # Bootstrap with 10000 resamples
    boot_result <- boot(trt_data, ratio_mean, R = 10000)

    # Calculate confidence intervals
    ci_result <- boot.ci(boot_result, type = "perc", conf = 0.95)

    bootstrap_results[[trt]] <- data.frame(
      Treatment = trt,
      Mean_Ratio = round(mean(trt_data$worms_per_infected, na.rm = TRUE), 3),
      Bootstrap_Mean = round(boot_result$t0, 3),
      CI_Lower = round(ci_result$percent[4], 3),
      CI_Upper = round(ci_result$percent[5], 3),
      N_Tanks = nrow(trt_data)
    )
  } else {
    # For treatments with only one observation, just report the mean
    bootstrap_results[[trt]] <- data.frame(
      Treatment = trt,
      Mean_Ratio = round(mean(trt_data$worms_per_infected, na.rm = TRUE), 3),
      Bootstrap_Mean = round(mean(trt_data$worms_per_infected, na.rm = TRUE), 3),
      CI_Lower = NA,
      CI_Upper = NA,
      N_Tanks = nrow(trt_data)
    )
  }
}

# Combine bootstrap results
bootstrap_summary <- do.call(rbind, bootstrap_results) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Bootstrap Confidence Intervals - Normalized Worm Burden",
    subtitle = "95% confidence intervals for worms per infected fish by treatment (tank-level data)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(Mean_Ratio, Bootstrap_Mean, CI_Lower, CI_Upper),
    decimals = 3
  )

# Method 5: Pairwise permutation tests
set.seed(42)
pairwise_permutation_results <- list()
treatments <- unique(tank_normalized_for_test$Treatment)

for(i in 1:(length(treatments)-1)) {
  for(j in (i+1):length(treatments)) {
    trt1 <- treatments[i]
    trt2 <- treatments[j]

    # Get data for these two treatments
    data1 <- tank_normalized_for_test %>% dplyr::filter(Treatment == trt1)
    data2 <- tank_normalized_for_test %>% dplyr::filter(Treatment == trt2)

    # Combine data for permutation test
    combined_data <- rbind(data1, data2)

    # Perform permutation test
    perm_test <- coin::oneway_test(
      worms_per_infected ~ Treatment,
      data = combined_data,
      distribution = approximate(nresample = 10000)
    )

    pairwise_permutation_results[[paste(trt1, "vs", trt2)]] <- data.frame(
      Comparison = paste(trt1, "vs", trt2),
      Z_statistic = round(statistic(perm_test), 4),
      P_value = round(pvalue(perm_test), 6),
      Significant = ifelse(pvalue(perm_test) < 0.05, "Yes", "No")
    )
  }
}

# Combine pairwise permutation results
pairwise_permutation_summary <- do.call(rbind, pairwise_permutation_results) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pairwise Permutation Tests - Normalized Worm Burden",
    subtitle = "Non-parametric pairwise comparisons for ratio data (tank-level data)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )

# Model comparison and diagnostics
normalized_model_comparison <- data.frame(
  Test = c("ANOVA", "Kruskal-Wallis", "Permutation Test"),
  Statistic = c(
    paste("F =", round(f_stat, 4)),
    paste("χ² =", round(normalized_kruskal$statistic, 4)),
    paste("Z =", round(statistic(normalized_permutation), 4))
  ),
  P_value = c(
    round(p_val, 6),
    round(normalized_kruskal$p.value, 6),
    round(pvalue(normalized_permutation), 6)
  ),
  Assumptions = c(
    "Normality, homogeneity of variance",
    "None (non-parametric)",
    "None (non-parametric)"
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Model Comparison - Normalized Worm Burden Analysis",
    subtitle = "Comparing different statistical approaches for ratio data"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Tank-level normalized worm burden analysis
set.seed(42)

# Tank-level normalized worm burden summary
tank_normalized_summary <- tank_normalized_data %>%
  dplyr::select(Treatment, Tank.ID, total_worms, n_infected, worms_per_infected, n_fish, infection_rate) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Tank-Level Normalized Worm Burden Summary",
    subtitle = "Worms per infected fish by treatment and tank"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(total_worms, n_infected, n_fish),
    decimals = 0
  ) %>%
  gt::fmt_number(
    columns = c(worms_per_infected, infection_rate),
    decimals = 2
  )

# Statistical test for tank-level normalized worm burden
# Using ANOVA on worms per infected fish (excluding tanks with 0 infected fish)
set.seed(42)

# One-way ANOVA for worms per infected fish across treatments
tank_normalized_anova <- aov(worms_per_infected ~ Treatment, data = tank_normalized_for_test)
tank_normalized_anova_summary <- summary(tank_normalized_anova)

# Extract ANOVA results
tank_normalized_anova_results <- data.frame(
  Statistic = c("F-value", "Degrees of freedom (Treatment)", "Degrees of freedom (Residuals)", "P-value"),
  Value = c(
    round(tank_normalized_anova_summary[[1]]$"F value"[1], 4),
    tank_normalized_anova_summary[[1]]$"Df"[1],
    tank_normalized_anova_summary[[1]]$"Df"[2],
    round(tank_normalized_anova_summary[[1]]$"Pr(>F)"[1], 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "ANOVA Results - Tank-Level Normalized Worm Burden",
    subtitle = "One-way ANOVA: worms per infected fish ~ Treatment"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Post-hoc Tukey test for tank-level normalized worm burden
set.seed(42)
tank_normalized_tukey <- TukeyHSD(tank_normalized_anova)

# Convert Tukey results to data frame
tank_normalized_tukey_df <- as.data.frame(tank_normalized_tukey$Treatment) %>%
  dplyr::mutate(
    comparison = rownames(.),
    diff = round(diff, 4),
    lwr = round(lwr, 4),
    upr = round(upr, 4),
    p_adj = round(`p adj`, 6),
    p_signif = ifelse(p_adj < 0.001, "***",
                     ifelse(p_adj < 0.01, "**",
                            ifelse(p_adj < 0.05, "*", "ns")))
  ) %>%
  dplyr::select(comparison, diff, lwr, upr, p_adj, p_signif)

# Display Tukey test results
tank_normalized_tukey_summary <- tank_normalized_tukey_df %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Tukey HSD Test Results - Tank-Level Normalized Worm Burden",
    subtitle = "Post-hoc pairwise comparisons for worms per infected fish"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(diff, lwr, upr, p_adj),
    decimals = 3
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = p_signif,
      rows = p_signif != "ns"
    )
  )

Table

Infection Prevalence

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Summary

Infection Prevalence Contingency Table
Number of infected vs uninfected fish by treatment
Treatment	n_infected	n_uninfected	n_fish	percent_infected	percent_uninfected
A- T- P+	31	51	82	37.8	62.2
A+ T- P+	19	61	80	23.8	76.2
A- T+ P+	21	56	77	27.3	72.7
A+ T+ P+	15	55	70	21.4	78.6

Infection Chi-Square

Chi-Square Test Results - Infection Prevalence
Test for independence between treatment and infection status
Statistic	Value
Chi-square statistic	6.16510
Degrees of freedom	3.00000
P-value	0.10385

Infection Pairwise

Pairwise Chi-Square Tests - Infection Prevalence
Post-hoc comparisons between treatment pairs
Comparison	Chi_square	P_value	Significant
A- T- P+ vs A+ T- P+	3.1190	0.077384	No
A- T- P+ vs A- T+ P+	1.5515	0.212910	No
A- T- P+ vs A+ T+ P+	4.0541	0.044064	Yes
A+ T- P+ vs A- T+ P+	0.1045	0.746535	No
A+ T- P+ vs A+ T+ P+	0.0205	0.886027	No
A- T+ P+ vs A+ T+ P+	0.3980	0.528098	No

Worm Burden

Summary

Worm Burden Summary by Treatment
Total, mean, median, standard deviation, and standard error of worm counts per treatment
Treatment	total_worms	mean_worms	median_worms	sd_worms	se_worms	n_fish
A- T- P+	167	2.04	0.00	5.16	0.57	82
A+ T- P+	45	0.56	0.00	1.58	0.18	80
A- T+ P+	121	1.57	0.00	3.59	0.41	77
A+ T+ P+	77	1.10	0.00	2.82	0.34	70

Negative Binomial Model

Negative Binomial Regression Results - Worm Burden
Model: glm.nb(worm counts ~ Treatment)
Treatment	Estimate	Std_Error	z_value	p_value	p_signif
Intercept	0.711	0.292	2.433	0.015	*
A+ T- P+	−1.287	0.435	−2.958	0.003	**
A- T+ P+	−0.259	0.422	−0.614	0.539	ns
A+ T+ P+	−0.616	0.438	−1.407	0.159	ns

Model Fit Stats

Negative Binomial Model Fit Statistics
Model diagnostics and fit measures
Statistic	Value
Log-likelihood	-388.0400
AIC	786.0800
Theta	0.1534
SE Theta	0.0226
Deviance	189.0100
Residual Deviance	189.0100

Worm Burden Pairwise

Pairwise Comparisons - Worm Burden (Negative Binomial)
Bonferroni-adjusted pairwise comparisons between treatments
contrast	estimate	SE	z.ratio	p.value	p_signif
(A- T- P+) - (A+ T- P+)	1.287	0.435	2.958	0.019	*
(A- T- P+) - (A- T+ P+)	0.259	0.422	0.614	1.000	ns
(A- T- P+) - (A+ T+ P+)	0.616	0.438	1.407	0.956	ns
(A+ T- P+) - (A- T+ P+)	−1.027	0.443	−2.317	0.123	ns
(A+ T- P+) - (A+ T+ P+)	−0.671	0.458	−1.464	0.859	ns
(A- T+ P+) - (A+ T+ P+)	0.357	0.446	0.800	1.000	ns

Model Diagnostics

Negative Binomial Model Diagnostics
Checking model adequacy and dispersion
Statistic	Value
Theta (dispersion parameter)	0.1534
SE Theta	0.0226
Deviance	189.01
Degrees of Freedom	305
Deviance/DF	0.62
Model adequacy	Good fit

Model Comparison

Model Comparison: Poisson vs Negative Binomial
Comparing model fit and adequacy
Model	Log_Likelihood	AIC	Deviance	DF	Deviance_DF_Ratio
Poisson	−800.23	1,608.45	1,338.72	305	4.39
Negative Binomial	−388.04	786.08	189.01	305	0.62

Normalized Worm Burden

Summary

Normalized Worm Burden Summary by Treatment
Worms per infected fish (normalized for infection prevalence)
Treatment	total_worms	n_infected	worms_per_infected	n_fish	infection_rate
A- T- P+	167	31	5.39	82	37.80
A+ T- P+	45	19	2.37	80	23.80
A- T+ P+	121	21	5.76	77	27.30
A+ T+ P+	77	15	5.13	70	21.40

ANOVA Results

ANOVA Results - Normalized Worm Burden
One-way ANOVA: worms per infected fish ~ Treatment (tank-level data)
Statistic	Value
F-value	1.857900
Degrees of freedom (Treatment)	3.000000
Degrees of freedom (Residuals)	8.000000
P-value	0.215039

Tukey HSD Test

Tukey HSD Test Results - Normalized Worm Burden
Post-hoc pairwise comparisons for worms per infected fish (tank-level data)
comparison	diff	lwr	upr	p_adj	p_signif
A+ T- P+-A- T- P+	−3.012	−8.471	2.448	0.354	ns
A- T+ P+-A- T- P+	0.198	−5.262	5.657	0.999	ns
A+ T+ P+-A- T- P+	0.536	−4.924	5.995	0.988	ns
A- T+ P+-A+ T- P+	3.210	−2.250	8.669	0.307	ns
A+ T+ P+-A+ T- P+	3.548	−1.912	9.007	0.238	ns
A+ T+ P+-A- T+ P+	0.338	−5.122	5.798	0.997	ns

Kruskal-Wallis Test

Kruskal-Wallis Test Results - Normalized Worm Burden
Non-parametric test for differences in worms per infected fish across treatments (tank-level data)
Statistic	Value
Kruskal-Wallis chi-squared	3.769200
Degrees of freedom	3.000000
P-value	0.287485

Permutation Test

Permutation Test Results - Normalized Worm Burden
Non-parametric permutation test for ratio data (tank-level data)
Statistic	Value
Z statistic	4.5168
P-value (approximate)	0.2356
Number of permutations	10000.0000

Bootstrap Confidence Intervals

Bootstrap Confidence Intervals - Normalized Worm Burden
95% confidence intervals for worms per infected fish by treatment (tank-level data)
Treatment	Mean_Ratio	Bootstrap_Mean	CI_Lower	CI_Upper	N_Tanks
A- T- P+	5.095	5.095	3.286	6.167	3
A+ T- P+	2.083	2.083	1.000	3.625	3
A- T+ P+	5.293	5.293	2.333	7.545	3
A+ T+ P+	5.631	5.631	3.143	8.000	3

Pairwise Permutation Tests

Pairwise Permutation Tests - Normalized Worm Burden
Non-parametric pairwise comparisons for ratio data (tank-level data)
Comparison	Z_statistic	P_value	Significant
A- T- P+ vs A+ T- P+	1.7453	0.2013	No
A- T- P+ vs A- T+ P+	-0.1230	0.9007	No
A- T- P+ vs A+ T+ P+	-0.3536	1.0000	No
A+ T- P+ vs A- T+ P+	-1.5176	0.2005	No
A+ T- P+ vs A+ T+ P+	-1.6551	0.2003	No
A- T+ P+ vs A+ T+ P+	-0.1804	0.8006	No

Model Comparison

Model Comparison - Normalized Worm Burden Analysis
Comparing different statistical approaches for ratio data
Test	Statistic	P_value	Assumptions
ANOVA	F = 1.8579	0.215039	Normality, homogeneity of variance
Kruskal-Wallis	χ² = 3.7692	0.287485	None (non-parametric)
Permutation Test	Z = 4.5168	0.235600	None (non-parametric)

Tank-Level Analysis

Tank-Level Normalized Worm Burden Summary
Worms per infected fish by treatment and tank
Treatment	Tank.ID	total_worms	n_infected	worms_per_infected	n_fish	infection_rate
A- T- P+	4	70	12	5.83	29	41.40
A- T- P+	5	23	7	3.29	25	28.00
A- T- P+	6	74	12	6.17	28	42.90
A+ T- P+	16	13	8	1.62	29	27.60
A+ T- P+	17	29	8	3.62	28	28.60
A+ T- P+	18	3	3	1.00	23	13.00
A- T+ P+	10	24	4	6.00	21	19.00
A- T+ P+	11	14	6	2.33	30	20.00
A- T+ P+	12	83	11	7.55	26	42.30
A+ T+ P+	22	22	7	3.14	28	25.00
A+ T+ P+	23	23	4	5.75	23	17.40
A+ T+ P+	24	32	4	8.00	19	21.10

ANOVA Results - Tank-Level Normalized Worm Burden
One-way ANOVA: worms per infected fish ~ Treatment
Statistic	Value
F-value	1.857900
Degrees of freedom (Treatment)	3.000000
Degrees of freedom (Residuals)	8.000000
P-value	0.215039

Tukey HSD Test Results - Tank-Level Normalized Worm Burden
Post-hoc pairwise comparisons for worms per infected fish
comparison	diff	lwr	upr	p_adj	p_signif
A+ T- P+-A- T- P+	−3.012	−8.471	2.448	0.354	ns
A- T+ P+-A- T- P+	0.198	−5.262	5.657	0.999	ns
A+ T+ P+-A- T- P+	0.536	−4.924	5.995	0.988	ns
A- T+ P+-A+ T- P+	3.210	−2.250	8.669	0.307	ns
A+ T+ P+-A+ T- P+	3.548	−1.912	9.007	0.238	ns
A+ T+ P+-A- T+ P+	0.338	−5.122	5.798	0.997	ns

03) Is host-microbiome recovery historically contingent?

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality

All

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

By Tank

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

Statistical Analysis - Recovery Mortality

Code

# Prepare data for statistical analysis - recovery mortality
recovery_mortality_data <- data.mortality %>%
  dplyr::filter(Treatment %in% c("A- T- P-", "A+ T- P-", "A- T+ P-", "A+ T+ P-")) %>%
  dplyr::group_by(Treatment) %>%
  dplyr::count() %>%
  dplyr::mutate(
    total_per_group = 90,
    dead = total_per_group - n,
    percent_mortality = round((dead / total_per_group) * 100, 1)
  )

# Chi-square test for independence between treatment and mortality
set.seed(42)
recovery_chi_square_test <- chisq.test(
  matrix(c(recovery_mortality_data$dead, recovery_mortality_data$n), 
         ncol = 2, byrow = FALSE)
)

# Create mortality rate summary table for display
recovery_mortality_rate_table <- recovery_mortality_data %>%
  dplyr::select(Treatment, dead, n, percent_mortality) %>%
  dplyr::rename(Survived = n) %>%
  dplyr::mutate(
    Total = dead + Survived,
    percent_survived = round((Survived / Total) * 100, 1)
  ) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Mortality Rate Summary - Recovery",
    subtitle = "Percent mortality and survival rates by treatment (parasite-unexposed only)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(dead, Survived, Total),
    decimals = 0
  ) %>%
  gt::fmt_number(
    columns = c(percent_mortality, percent_survived),
    decimals = 1
  )

# Chi-square test results
recovery_chi_square_results <- data.frame(
  Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
  Value = c(
    round(recovery_chi_square_test$statistic, 4),
    recovery_chi_square_test$parameter,
    round(recovery_chi_square_test$p.value, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Chi-Square Test Results - Recovery Mortality",
    subtitle = "Test for independence between treatment and mortality (parasite-unexposed only)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Post-hoc pairwise comparisons using chi-square tests
set.seed(42)
recovery_pairwise_results <- list()
recovery_treatments <- recovery_mortality_data$Treatment

for(i in 1:(length(recovery_treatments)-1)) {
  for(j in (i+1):length(recovery_treatments)) {
    trt1 <- recovery_treatments[i]
    trt2 <- recovery_treatments[j]
    
    # Get data for these two treatments
    data1 <- recovery_mortality_data %>% dplyr::filter(Treatment == trt1)
    data2 <- recovery_mortality_data %>% dplyr::filter(Treatment == trt2)
    
    # Create 2x2 contingency table
    cont_table <- matrix(
      c(data1$dead, data1$n, data2$dead, data2$n),
      nrow = 2, ncol = 2,
      dimnames = list(
        c("Dead", "Survived"),
        c(trt1, trt2)
      )
    )
    
    # Perform chi-square test
    test_result <- chisq.test(cont_table)
    
    recovery_pairwise_results[[paste(trt1, "vs", trt2)]] <- data.frame(
      Comparison = paste(trt1, "vs", trt2),
      Chi_square = round(test_result$statistic, 4),
      P_value = round(test_result$p.value, 6),
      Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
    )
  }
}

# Combine pairwise results
recovery_pairwise_summary <- do.call(rbind, recovery_pairwise_results) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pairwise Chi-Square Tests - Recovery Mortality",
    subtitle = "Post-hoc comparisons between treatment pairs (parasite-unexposed only)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )

Table

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality Rates

Mortality Rate Summary - Recovery
Percent mortality and survival rates by treatment (parasite-unexposed only)
dead	Survived	percent_mortality	Total	percent_survived
A- T- P-
9	81	10.0	90	90.0
A+ T- P-
14	76	15.6	90	84.4
A- T+ P-
19	71	21.1	90	78.9
A+ T+ P-
30	60	33.3	90	66.7

Chi-Square

Chi-Square Test Results - Recovery Mortality
Test for independence between treatment and mortality (parasite-unexposed only)
Statistic	Value
Chi-square statistic	16.805600
Degrees of freedom	3.000000
P-value	0.000775

Chi-Square Pairwise

Pairwise Chi-Square Tests - Recovery Mortality
Post-hoc comparisons between treatment pairs (parasite-unexposed only)
Comparison	Chi_square	P_value	Significant
A- T- P- vs A+ T- P-	0.7976	0.371823	No
A- T- P- vs A- T+ P-	3.4258	0.064187	No
A- T- P- vs A+ T+ P-	13.0933	0.000296	Yes
A+ T- P- vs A- T+ P-	0.5937	0.440995	No
A+ T- P- vs A+ T+ P-	6.7680	0.009280	Yes
A- T+ P- vs A+ T+ P-	2.8042	0.094019	No

04) History of stressors

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality

All

By Tank

Statistical Analysis - History of Stressors Mortality

Code

# Prepare data for statistical analysis - percent mortality rates by History
history_mortality_stats_data <- data.mortality %>%
  dplyr::group_by(History) %>%
  dplyr::count() %>%
  dplyr::mutate(
    # Calculate expected totals: 90 fish per treatment group
    total_per_group = dplyr::case_when(
      History == 0 ~ 90 * 2,  # 2 treatment groups: A- T- P- and A- T- P+
      History == 1 ~ 90 * 4,  # 4 treatment groups: A+ T- P-, A+ T- P+, A- T+ P-, A- T+ P+
      History == 2 ~ 90 * 2,  # 2 treatment groups: A+ T+ P- and A+ T+ P+
      TRUE ~ 0
    ),
    dead = total_per_group - n,
    percent_mortality = round((dead / total_per_group) * 100, 1),
    History_Label = dplyr::case_when(
      History == 0 ~ "No prior stressors",
      History == 1 ~ "One prior stressor", 
      History == 2 ~ "Two prior stressors",
      TRUE ~ "Unknown"
    )
  )

# Chi-square test for independence between history and mortality
set.seed(42)
history_chi_square_test <- chisq.test(
  matrix(c(history_mortality_stats_data$dead, history_mortality_stats_data$n), 
         ncol = 2, byrow = FALSE)
)

# Create mortality rate summary table for display
history_mortality_rate_table <- history_mortality_stats_data %>%
  dplyr::select(History, History_Label, dead, n, percent_mortality) %>%
  dplyr::rename(Survived = n) %>%
  dplyr::mutate(
    Total = dead + Survived,
    percent_survived = round((Survived / Total) * 100, 1)
  ) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Mortality Rate Summary - History of Stressors",
    subtitle = "Percent mortality and survival rates by history of stressors"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(dead, Survived, Total),
    decimals = 0
  ) %>%
  gt::fmt_number(
    columns = c(percent_mortality, percent_survived),
    decimals = 1
  )

# Chi-square test results
history_test_stat_table <- data.frame(
  Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
  Value = c(
    round(history_chi_square_test$statistic, 4),
    history_chi_square_test$parameter,
    round(history_chi_square_test$p.value, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Chi-Square Test Results - History of Stressors Mortality",
    subtitle = "Test for independence between history of stressors and mortality"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Post-hoc pairwise comparisons using chi-square tests
set.seed(42)
history_pairwise_results <- list()
history_levels <- history_mortality_stats_data$History

for(i in 1:(length(history_levels)-1)) {
  for(j in (i+1):length(history_levels)) {
    hist1 <- history_levels[i]
    hist2 <- history_levels[j]
    
    # Get data for these two history levels
    data1 <- history_mortality_stats_data %>% dplyr::filter(History == hist1)
    data2 <- history_mortality_stats_data %>% dplyr::filter(History == hist2)
    
    # Create 2x2 contingency table
    cont_table <- matrix(
      c(data1$dead, data1$n, data2$dead, data2$n),
      nrow = 2, ncol = 2,
      dimnames = list(
        c("Dead", "Survived"),
        c(paste("History", hist1), paste("History", hist2))
      )
    )
    
    # Perform chi-square test
    test_result <- chisq.test(cont_table)
    
    history_pairwise_results[[paste("History", hist1, "vs", "History", hist2)]] <- data.frame(
      Comparison = paste("History", hist1, "vs", "History", hist2),
      Chi_square = round(test_result$statistic, 4),
      P_value = round(test_result$p.value, 6),
      Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
    )
  }
}

# Combine pairwise results
history_pairwise_summary <- do.call(rbind, history_pairwise_results) %>%
  dplyr::arrange(desc(Significant), P_value) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pairwise Chi-Square Tests - History of Stressors Mortality",
    subtitle = "Post-hoc comparisons between history levels"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )

# Show only significant pairwise comparisons
history_pairwise_summary__SigOnly <- do.call(rbind, history_pairwise_results) %>%
  dplyr::arrange(desc(Significant), P_value) %>%
  dplyr::filter(Significant == "Yes") %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Significant Pairwise Chi-Square Tests - History of Stressors Mortality",
    subtitle = "Post-hoc comparisons between history levels (significant only)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )

# Additional analysis: Trend analysis across history levels
# Calculate trend in mortality rates
history_trend_data <- history_mortality_stats_data %>%
  dplyr::arrange(History) %>%
  dplyr::select(History, percent_mortality)

# Linear trend test (Cochran-Armitage trend test equivalent using chi-square)
set.seed(42)
# Create trend weights (0, 1, 2 for history levels)
trend_weights <- c(0, 1, 2)
trend_chi_square <- chisq.test(
  matrix(c(history_mortality_stats_data$dead, history_mortality_stats_data$n), 
         ncol = 2, byrow = FALSE),
  simulate.p.value = TRUE
)

# Calculate correlation between history level and mortality rate
history_correlation <- cor.test(
  history_trend_data$History, 
  history_trend_data$percent_mortality, 
  method = "pearson"
)

# Trend analysis results
history_trend_results <- data.frame(
  Statistic = c("Correlation coefficient (r)", "P-value (correlation)", "Chi-square trend test p-value"),
  Value = c(
    round(history_correlation$estimate, 4),
    round(history_correlation$p.value, 6),
    round(trend_chi_square$p.value, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Trend Analysis - History of Stressors Mortality",
    subtitle = "Testing for linear trend in mortality rates across history levels"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Effect size analysis: Calculate odds ratios between history levels
history_effect_sizes <- data.frame(
  Comparison = c("History 1 vs History 0", "History 2 vs History 0", "History 2 vs History 1"),
  Odds_Ratio = c(
    round((history_mortality_stats_data$dead[2] * history_mortality_stats_data$n[1]) / 
          (history_mortality_stats_data$dead[1] * history_mortality_stats_data$n[2]), 3),
    round((history_mortality_stats_data$dead[3] * history_mortality_stats_data$n[1]) / 
          (history_mortality_stats_data$dead[1] * history_mortality_stats_data$n[3]), 3),
    round((history_mortality_stats_data$dead[3] * history_mortality_stats_data$n[2]) / 
          (history_mortality_stats_data$dead[2] * history_mortality_stats_data$n[3]), 3)
  ),
  Mortality_Rate_1 = c(
    round(history_mortality_stats_data$percent_mortality[2], 1),
    round(history_mortality_stats_data$percent_mortality[3], 1),
    round(history_mortality_stats_data$percent_mortality[3], 1)
  ),
  Mortality_Rate_2 = c(
    round(history_mortality_stats_data$percent_mortality[1], 1),
    round(history_mortality_stats_data$percent_mortality[1], 1),
    round(history_mortality_stats_data$percent_mortality[2], 1)
  )
) %>%
  dplyr::mutate(
    Interpretation = dplyr::case_when(
      Odds_Ratio > 1 ~ "Higher mortality in first group",
      Odds_Ratio < 1 ~ "Lower mortality in first group",
      Odds_Ratio == 1 ~ "No difference"
    )
  ) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Effect Size Analysis - History of Stressors Mortality",
    subtitle = "Odds ratios comparing mortality rates between history levels"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(Odds_Ratio, Mortality_Rate_1, Mortality_Rate_2),
    decimals = 1
  )

Table

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality Rates

Mortality Rate Summary - History of Stressors
Percent mortality and survival rates by history of stressors
History_Label	dead	Survived	percent_mortality	Total	percent_survived
0
No prior stressors	17	163	9.4	180	90.6
1
One prior stressor	56	304	15.6	360	84.4
2
Two prior stressors	50	130	27.8	180	72.2

Chi-Square

Chi-Square Test Results - History of Stressors Mortality
Test for independence between history of stressors and mortality
Statistic	Value
Chi-square statistic	22.542000
Degrees of freedom	2.000000
P-value	0.000013

Chi-Square Pairwise

Pairwise Chi-Square Tests - History of Stressors Mortality
Post-hoc comparisons between history levels
Comparison	Chi_square	P_value	Significant
History 0 vs History 2	18.7785	0.000015	Yes
History 1 vs History 2	10.6010	0.001130	Yes
History 0 vs History 1	3.3284	0.068094	No

Significant Pairwise Chi-Square Tests - History of Stressors Mortality
Post-hoc comparisons between history levels (significant only)
Comparison	Chi_square	P_value	Significant
History 0 vs History 2	18.7785	0.000015	Yes
History 1 vs History 2	10.6010	0.001130	Yes

Trend Analysis

Trend Analysis - History of Stressors Mortality
Testing for linear trend in mortality rates across history levels
Statistic	Value
Correlation coefficient (r)	0.982700
P-value (correlation)	0.118467
Chi-square trend test p-value	0.000500

Effect Sizes

Effect Size Analysis - History of Stressors Mortality
Odds ratios comparing mortality rates between history levels
Comparison	Odds_Ratio	Mortality_Rate_1	Mortality_Rate_2	Interpretation
History 1 vs History 0	1.8	15.6	9.4	Higher mortality in first group
History 2 vs History 0	3.7	27.8	9.4	Higher mortality in first group
History 2 vs History 1	2.1	27.8	15.6	Higher mortality in first group

05) History of stressors by pathogen exposure

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality

All

By Tank

Statistical Analysis - History × Pathogen Interaction

Code

# Prepare data for statistical analysis - History × Pathogen interaction
interaction_stats_data <- history_pathogen_interaction_data %>%
  dplyr::select(History, Pathogen, History_Label, Pathogen_Label, dead, n, percent_mortality) %>%
  dplyr::rename(Survived = n) %>%
  dplyr::mutate(
    Total = dead + Survived,
    percent_survived = round((Survived / Total) * 100, 1)
  )

# 1. Chi-square test for independence between History and Pathogen exposure
set.seed(42)
# Create 3x2 contingency table (History × Pathogen)
interaction_contingency <- matrix(
  c(interaction_stats_data$dead), 
  nrow = 3, ncol = 2,
  dimnames = list(
    c("History 0", "History 1", "History 2"),
    c("Unexposed", "Exposed")
  )
)

interaction_chi_square <- chisq.test(interaction_contingency)

# 2. Test for interaction effect using logistic regression approach
# Create individual-level data for interaction analysis
set.seed(42)
individual_interaction_data <- data.mortality %>%
  dplyr::select(History, Pathogen, Sample) %>%
  dplyr::mutate(
    # Create binary outcome: 1 = survived, 0 = died (assuming all sampled fish survived)
    Outcome = 1,  # All fish in the dataset survived (were sampled)
    History = factor(History, levels = c(0, 1, 2)),
    Pathogen = factor(Pathogen, levels = c(0, 1))
  )

# Fit logistic regression model with interaction
interaction_model <- glm(Outcome ~ History * Pathogen, 
                        data = individual_interaction_data, 
                        family = "binomial")

## Warning: glm.fit: algorithm did not converge

# Get model summary
interaction_summary <- summary(interaction_model)

# Extract coefficients and significance
interaction_coefficients <- data.frame(
  Term = rownames(interaction_summary$coefficients),
  Estimate = round(interaction_summary$coefficients[, 1], 4),
  Std_Error = round(interaction_summary$coefficients[, 2], 4),
  z_value = round(interaction_summary$coefficients[, 3], 4),
  p_value = round(interaction_summary$coefficients[, 4], 6)
) %>%
  dplyr::mutate(
    p_signif = ifelse(p_value < 0.001, "***", 
                     ifelse(p_value < 0.01, "**",
                            ifelse(p_value < 0.05, "*", "ns")))
  )

# 3. Create summary tables for display
interaction_summary_table <- interaction_stats_data %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Mortality Summary - History × Pathogen Interaction",
    subtitle = "Percent mortality and survival rates by history and pathogen exposure"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(dead, Survived, Total),
    decimals = 0
  ) %>%
  gt::fmt_number(
    columns = c(percent_mortality, percent_survived),
    decimals = 1
  )

# Chi-square test results
interaction_chi_results <- data.frame(
  Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
  Value = c(
    round(interaction_chi_square$statistic, 4),
    interaction_chi_square$parameter,
    round(interaction_chi_square$p.value, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Chi-Square Test Results - History × Pathogen Interaction",
    subtitle = "Test for independence between history and pathogen exposure"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Logistic regression results
interaction_model_results <- interaction_coefficients %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Logistic Regression Results - History × Pathogen Interaction",
    subtitle = "Model: Outcome ~ History * Pathogen"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(Estimate, Std_Error, z_value, p_value),
    decimals = 3
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = p_signif,
      rows = p_signif != "ns"
    )
  )

# 4. Calculate effect sizes for each comparison
# Calculate odds ratios for each history level comparing exposed vs unexposed
set.seed(42)
effect_size_comparisons <- data.frame(
  Comparison = c(
    "History 0: Exposed vs Unexposed",
    "History 1: Exposed vs Unexposed", 
    "History 2: Exposed vs Unexposed"
  ),
  Odds_Ratio = c(
    # History 0: Exposed vs Unexposed
    round((interaction_stats_data$dead[2] * interaction_stats_data$Survived[1]) / 
          (interaction_stats_data$dead[1] * interaction_stats_data$Survived[2]), 3),
    # History 1: Exposed vs Unexposed
    round((interaction_stats_data$dead[4] * interaction_stats_data$Survived[3]) / 
          (interaction_stats_data$dead[3] * interaction_stats_data$Survived[4]), 3),
    # History 2: Exposed vs Unexposed
    round((interaction_stats_data$dead[6] * interaction_stats_data$Survived[5]) / 
          (interaction_stats_data$dead[5] * interaction_stats_data$Survived[6]), 3)
  ),
  Mortality_Exposed = c(
    round(interaction_stats_data$percent_mortality[2], 1),
    round(interaction_stats_data$percent_mortality[4], 1),
    round(interaction_stats_data$percent_mortality[6], 1)
  ),
  Mortality_Unexposed = c(
    round(interaction_stats_data$percent_mortality[1], 1),
    round(interaction_stats_data$percent_mortality[3], 1),
    round(interaction_stats_data$percent_mortality[5], 1)
  )
) %>%
  dplyr::mutate(
    Mortality_Difference = Mortality_Exposed - Mortality_Unexposed,
    Interpretation = dplyr::case_when(
      Odds_Ratio > 1 ~ "Higher mortality when exposed",
      Odds_Ratio < 1 ~ "Lower mortality when exposed",
      Odds_Ratio == 1 ~ "No difference"
    )
  ) 

effect_size_comparisons_table <- effect_size_comparisons %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Effect Size Analysis - History × Pathogen Interaction",
    subtitle = "Odds ratios comparing exposed vs unexposed within each history level"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(Odds_Ratio, Mortality_Exposed, Mortality_Unexposed, Mortality_Difference),
    decimals = 1
  )

# 5. Test for trend in pathogen effect across history levels
set.seed(42)
# Calculate correlation between history level and the difference in mortality (exposed - unexposed)
trend_data <- effect_size_comparisons %>%
  dplyr::mutate(History_Level = c(0, 1, 2))

pathogen_effect_trend <- cor.test(
  trend_data$History_Level, 
  trend_data$Mortality_Difference, 
  method = "pearson"
)

# Trend analysis results
pathogen_effect_trend_results <- data.frame(
  Statistic = c("Correlation coefficient (r)", "P-value", "Interpretation"),
  Value = c(
    round(pathogen_effect_trend$estimate, 4),
    round(pathogen_effect_trend$p.value, 6),
    ifelse(pathogen_effect_trend$p.value < 0.05, 
           "Significant trend in pathogen effect across history levels",
           "No significant trend in pathogen effect across history levels")
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Trend Analysis - Pathogen Effect Across History Levels",
    subtitle = "Testing whether the effect of pathogen exposure changes across history levels"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# 6. Post-hoc pairwise comparisons within each pathogen group
set.seed(42)
# Within unexposed group
unexposed_data <- interaction_stats_data %>% dplyr::filter(Pathogen == 0)
unexposed_pairwise <- list()

for(i in 1:(nrow(unexposed_data)-1)) {
  for(j in (i+1):nrow(unexposed_data)) {
    hist1 <- unexposed_data$History[i]
    hist2 <- unexposed_data$History[j]
    
    # Create 2x2 contingency table
    cont_table <- matrix(
      c(unexposed_data$dead[i], unexposed_data$Survived[i], 
        unexposed_data$dead[j], unexposed_data$Survived[j]),
      nrow = 2, ncol = 2,
      dimnames = list(
        c("Dead", "Survived"),
        c(paste("History", hist1), paste("History", hist2))
      )
    )
    
    # Perform chi-square test
    test_result <- chisq.test(cont_table)
    
    unexposed_pairwise[[paste("History", hist1, "vs", "History", hist2)]] <- data.frame(
      Comparison = paste("History", hist1, "vs", "History", hist2, "(Unexposed)"),
      Chi_square = round(test_result$statistic, 4),
      P_value = round(test_result$p.value, 6),
      Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
    )
  }
}

# Within exposed group
exposed_data <- interaction_stats_data %>% dplyr::filter(Pathogen == 1)
exposed_pairwise <- list()

for(i in 1:(nrow(exposed_data)-1)) {
  for(j in (i+1):nrow(exposed_data)) {
    hist1 <- exposed_data$History[i]
    hist2 <- exposed_data$History[j]
    
    # Create 2x2 contingency table
    cont_table <- matrix(
      c(exposed_data$dead[i], exposed_data$Survived[i], 
        exposed_data$dead[j], exposed_data$Survived[j]),
      nrow = 2, ncol = 2,
      dimnames = list(
        c("Dead", "Survived"),
        c(paste("History", hist1), paste("History", hist2))
      )
    )
    
    # Perform chi-square test
    test_result <- chisq.test(cont_table)
    
    exposed_pairwise[[paste("History", hist1, "vs", "History", hist2)]] <- data.frame(
      Comparison = paste("History", hist1, "vs", "History", hist2, "(Exposed)"),
      Chi_square = round(test_result$statistic, 4),
      P_value = round(test_result$p.value, 6),
      Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
    )
  }
}

# 7. Pairwise comparisons between same history types across pathogen exposures
set.seed(42)
same_history_pathogen_pairwise <- list()

# For each history level, compare unexposed vs exposed
for(hist_level in c(0, 1, 2)) {
  # Get data for this history level
  unexposed_row <- interaction_stats_data %>% dplyr::filter(History == hist_level & Pathogen == 0)
  exposed_row <- interaction_stats_data %>% dplyr::filter(History == hist_level & Pathogen == 1)
  
  if(nrow(unexposed_row) > 0 && nrow(exposed_row) > 0) {
    # Calculate percent mortality for each group
    unexposed_mortality <- unexposed_row$percent_mortality
    exposed_mortality <- exposed_row$percent_mortality
    
    # Calculate the difference in percent mortality
    mortality_difference <- exposed_mortality - unexposed_mortality
    
    # Perform t-test on percent mortality rates
    # Since we have single values, we'll use a simple comparison
    # For small sample sizes, we can use a z-test for proportions
    # Calculate standard error of the difference
    unexposed_se <- sqrt((unexposed_mortality * (100 - unexposed_mortality)) / unexposed_row$Total)
    exposed_se <- sqrt((exposed_mortality * (100 - exposed_mortality)) / exposed_row$Total)
    se_diff <- sqrt(unexposed_se^2 + exposed_se^2)
    
    # Calculate z-statistic
    z_stat <- mortality_difference / se_diff
    p_value <- 2 * (1 - pnorm(abs(z_stat)))
    
    same_history_pathogen_pairwise[[paste("History", hist_level)]] <- data.frame(
      Comparison = paste("History", hist_level, "(Unexposed) vs History", hist_level, "(Exposed)"),
      Unexposed_Mortality = round(unexposed_mortality, 1),
      Exposed_Mortality = round(exposed_mortality, 1),
      Mortality_Difference = round(mortality_difference, 1),
      Z_statistic = round(z_stat, 4),
      P_value = round(p_value, 6),
      Significant = ifelse(p_value < 0.05, "Yes", "No")
    )
  }
}

# Create separate tables for different types of pairwise comparisons

# Table 1: Within unexposed group comparisons
unexposed_pairwise_results <- do.call(rbind, unexposed_pairwise) %>%
  dplyr::arrange(desc(Significant), P_value) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pairwise Chi-Square Tests - Within Unexposed Group",
    subtitle = "Comparing history levels within unexposed fish"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )

# Table 2: Within exposed group comparisons
exposed_pairwise_results <- do.call(rbind, exposed_pairwise) %>%
  dplyr::arrange(desc(Significant), P_value) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pairwise Chi-Square Tests - Within Exposed Group",
    subtitle = "Comparing history levels within exposed fish"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )

# Create separate table for same history comparisons across pathogen exposures
same_history_pathogen_results <- do.call(rbind, same_history_pathogen_pairwise) %>%
  dplyr::arrange(desc(Significant), P_value) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pathogen Effect Within Each History Level",
    subtitle = "Comparing percent mortality between unexposed vs exposed within the same history of stressors"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(Unexposed_Mortality, Exposed_Mortality, Mortality_Difference),
    decimals = 1
  ) %>%
  gt::fmt_number(
    columns = c(Z_statistic, P_value),
    decimals = 4
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )

Table

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Summary Table

Mortality Summary - History × Pathogen Interaction
Percent mortality and survival rates by history and pathogen exposure
History_Label	Pathogen_Label	dead	Survived	percent_mortality	Total	percent_survived
0 - 0
No prior stressors	Unexposed	9	81	10.0	90	90.0
0 - 1
No prior stressors	Exposed	8	82	8.9	90	91.1
1 - 0
One prior stressor	Unexposed	33	147	18.3	180	81.7
1 - 1
One prior stressor	Exposed	23	157	12.8	180	87.2
2 - 0
Two prior stressors	Unexposed	30	60	33.3	90	66.7
2 - 1
Two prior stressors	Exposed	20	70	22.2	90	77.8

Chi-Square Test

Chi-Square Test Results - History × Pathogen Interaction
Test for independence between history and pathogen exposure
Statistic	Value
Chi-square statistic	18.392800
Degrees of freedom	2.000000
P-value	0.000101

Logistic Regression

Logistic Regression Results - History × Pathogen Interaction
Model: Outcome ~ History * Pathogen
Term	Estimate	Std_Error	z_value	p_value	p_signif
(Intercept)	26.566	39,569.344	0.001	0.999	ns
History1	0.000	49,279.608	0.000	1.000	ns
History2	0.000	60,658.574	0.000	1.000	ns
Pathogen1	0.000	55,788.569	0.000	1.000	ns
History1:Pathogen1	0.000	69,158.557	0.000	1.000	ns
History2:Pathogen1	0.000	83,891.969	0.000	1.000	ns

Effect Sizes

Effect Size Analysis - History × Pathogen Interaction
Odds ratios comparing exposed vs unexposed within each history level
Comparison	Odds_Ratio	Mortality_Exposed	Mortality_Unexposed	Mortality_Difference	Interpretation
History 0: Exposed vs Unexposed	0.9	8.9	10.0	−1.1	Lower mortality when exposed
History 1: Exposed vs Unexposed	0.7	12.8	18.3	−5.5	Lower mortality when exposed
History 2: Exposed vs Unexposed	0.6	22.2	33.3	−11.1	Lower mortality when exposed

Trend Analysis

Trend Analysis - Pathogen Effect Across History Levels
Testing whether the effect of pathogen exposure changes across history levels
Statistic	Value
Correlation coefficient (r)	-0.9976
P-value	0.044036
Interpretation	Significant trend in pathogen effect across history levels

Within Unexposed Group

Pairwise Chi-Square Tests - Within Unexposed Group
Comparing history levels within unexposed fish
Comparison	Chi_square	P_value	Significant
History 0 vs History 2 (Unexposed)	13.0933	0.000296	Yes
History 1 vs History 2 (Unexposed)	6.7314	0.009473	Yes
History 0 vs History 1 (Unexposed)	2.5693	0.108955	No

Within Exposed Group

Pairwise Chi-Square Tests - Within Exposed Group
Comparing history levels within exposed fish
Comparison	Chi_square	P_value	Significant
History 0 vs History 2 (Exposed)	5.1175	0.023686	Yes
History 1 vs History 2 (Exposed)	3.3228	0.068326	No
History 0 vs History 1 (Exposed)	0.5512	0.457833	No

Same History Across Pathogen Exposure

Pathogen Effect Within Each History Level
Comparing percent mortality between unexposed vs exposed within the same history of stressors
Comparison	Unexposed_Mortality	Exposed_Mortality	Mortality_Difference	Z_statistic	P_value	Significant
History 2 (Unexposed) vs History 2 (Exposed)	33.3	22.2	−11.1	−1.6759	0.0938	No
History 1 (Unexposed) vs History 1 (Exposed)	18.3	12.8	−5.5	−1.4440	0.1487	No
History 0 (Unexposed) vs History 0 (Exposed)	10.0	8.9	−1.1	−0.2523	0.8008	No